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Comparative study of high-order 
positivity-preserving WENO schemes 

By D.V. Kotov, H.C. Yee and B. Sjogreen 


1. Motivation and objectives 

In gas dynamics and magnetohydrodynamics flows, physically, the density p and the 
pressure p should both be positive. In a standard conservative numerical scheme, however, 
the computed internal energy is obtained by subtracting the kinetic energy from the total 
energy, resulting in a computed p that may be negative. Examples are problems in which 
the dominant energy is kinetic. Negative p may often emerge in computing blast waves. 
In such situations the computed eigenvalues of the Jacobian will become imaginary. 
Consequently, the initial value problem for the linearized system will be ill posed. This 
explains why failure of preserving positivity of density or pressure may cause blow-ups 
of the numerical algorithm. The adhoc methods in numerical strategy which modify 
the computed negative density and/or the computed negative pressure to be positive 
are neither a conservative cure nor a stable solution. Conservative positivity-preserving 
schemes are more appropriate for such flow problems. 

The ideas of Zhang & Shu (2012) and Hu et al. (2012) precisely address the afore- 
mentioned issue. Zhang & Shu constructed a new conservative positivity-preserving pro- 
cedure to preserve positive density and pressure for high-order WENO schemes by the 
Lax- Friedrichs flux (WENO/LLF). In general, WENO/LLF is too dissipative for flows 
such as turbulence with strong shocks computed in direct numerical simulations (DNS) 
and large eddy simulations (LES). The new conservative positivity-preserving procedure 
proposed in Hu et al. (2012) can be used with any high-order shock-capturing scheme, 
including high-order WENO schemes using the Roe’s flux (WENO/Roe). 

The goal of this study is to compare the results obtained by non-positivity-preserving 
methods with the recently developed positivity-preserving schemes for representative 
test cases. In particular the more difficult 3D Noh and Sedov problems are considered. 
These test cases are chosen because of the negative pressure/density most often exhib- 
ited by standard high-order shock- capturing schemes. The simulation of a hypersonic 
nonequilibrium viscous shock tube that is related to the NASA Electric Arc Shock Tube 
(EAST) is also included. EAST is a high-temperature and high Mach number viscous non- 
equilibrium flow consisting of 13 species. In addition, as most common shock-capturing 
schemes have been developed for problems without source terms, when applied to prob- 
lems with nonlinear and/or stiff source terms these methods can result in spurious solu- 
tions, even when solving a conservative system of equations with a conservative scheme. 
This kind of behavior can be observed even for a scalar case (LeVeque & Yee 1990) 
as well as for the case consisting of two species and one reaction (Wang et al. 2012). 
For further information concerning this issue see (LeVeque & Yee 1990; Griffiths et al. 
1992; Lafon & Yee 1996; Yee et al. 2012). This EAST example indicated that standard 
high-order shock-capturing methods exhibit instability of density /pressure in addition 
to grid-dependent discontinuity locations with insufficient grid points. The evaluation 
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of these test cases is based on the stability of the numerical schemes together with the 
accuracy of the obtained solutions. 


2. Positivity-preserving algorithms 

Here we briefly describe the positivity-preserving method of Hu et al. ( 2012 ). Readers 
are referred to Zhang & Shu ( 2012 ) for their positivity-preserving WENO schemes that 
are valid only for the Lax-Friedrichs flux formulation. Consider the Euler equations: 
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where p is the density, u is the velocity, m is the momentum, E is the total energy, p 
is the pressure, e is the internal energy, and 7 > 1 is a constant. The speed of sound is 
given by c = \fjp~j~p and the three eigenvalues of the Jacobian f'(w) are u — c, u and 
u + c. 

A general explicit k th - order conservative scheme with Euler-forward time integration 
for Eq. (2.1) can be written as 
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where the superscripts n and n + 1 represent the old and new time-steps, respectively, 
and A = At/ Ax, where At is the time-step size and Ax is the grid-step size. 

Eq. (2.3) can be rewritten as follows: 
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The positivity-preserving procedure of Hu et al. involves the first-order Lax-Friedrichs 
scheme with numerical flux: 
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Let wf F,± = wf =p 2Af^ . The positive density is first enforced by the following cut-off 


flux limiter for positive density: 

1. For all i initialize 0f +1 / 2 = 1,0 
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2 . If p (w+) < e p , solve 6>+_ 1/2 from (1 - 0+ +1/2 )p(wf F ’ + ) + d+ +1/2 p(w+) = e p . 

3. If p (w“ ) < e p , solve 0~ +1/2 from (1 - 0~ +1/2 )p( wf F ~) + 0~ +1/2 p( w“ ) = e p . 

4. Set 0 pp+ 1/2 = min(0f +1 / 2 , 0 i+1 ^ 2 ), f i+ i = (1 — ^ pp+ i/ 2 )^ i+ 1 + ® p ,j+ 1 / 2 ^+j‘ 

Here, e p = min{10 -13 , p^ in }, where is the minimum density in the initial condi- 
tion, f* +i is the limited flux, and 0 < 07^2 < 1 are the limiting factors corresponding 

to the two neighboring cells, which share the same flux f i+ i- After applying this flux 
limiter, Eq. (2.4) becomes 
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The positive pressure is further enforced by the following cut-off flux limiter for positive 
pressure: 


1. For all i initialize Of +1 j 2 = 1 , 0, 
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Again, e p = min{10 13 ,P^u n }, where is the minimum pressure in the initial 

condition, and f*_j\ is the further limited flux. With these limited fluxes, the original 


scheme (2.3) is modified as 


w" +1 


= w " — a (t*+i — t*-i) 


(2.7) 


The above new procedure can be applied at each sub-stage of a TVD Runge-Kutta 
(Shu & Osher 1988) method, which is a convex combination of Euler-forward time-steps. 
Note that any high-order numerical flux can be used in this formulation. 


3. Numerical results 

The numerical experiments include ID and 3D Noh problems, 3D Sedov blast, a Mach 
2000 jet test case, and the ID EAST hypersonic viscous nonequilibrium shock tube. The 
comparative study includes the following numerical schemes: 

• UPWIND - first-order upwind scheme using Roe average state; 

• TVD - second-order TVD (Yee 1989; Yee et ai 1990); 

• Standard fifth- and seventh-order WENO (WEN05 and WEN07) with local Lax- 
Friedrichs flux (WEN05-LLF and WEN07-LLF) and Roe flux (WEN05-Roe), see Jiang 
& Shu (1996); 

• WEN05fi+split and WEN07fi+split - high-order nonlinear filter counterparts of 
WEN05 and WEN07 using the wavelet flow sensor in conjunction with the Ducros et 
al. splitting of the inviscid flux derivatives (see Yee & Sjogreen (2007, 2010)); 

• Zhang & Shu positivity-preserving WEN05 and WEN07 local Lax-Friedrichs flux 
(WEN05P, WEN07P) and of fifth order with global Lax-Friedrichs flux (WENOGP); 

• Hu et al. positivity-preserving WEN05 with local Lax-Friedrichs flux (WEN05PH) 
and Roe’s flux (WEN05PH-Roe). Similarly for WEN05Pfi+split. 

Note that for the TVD scheme, Roe’s average state is employed. 

3.1. The Noh problem 

The first test case is the well-known ID and 3D spherical Noh implosion problem (Noh 
1987). The initial conditions are p m 1, p = 0 and u = unit vector directed toward the 
origin with 7 s= 5/3. In this problem an infinite-strength shock expands outward from 
the origin at a constant velocity of 1/3. The goal is to test the ability of the scheme to 
preserve spherical symmetry and produce the correct entropy jump for adiabatic shock 
compression. 

The results obtained using different schemes for the ID case are shown in Figure 1. The 
grid size is h — 0.002. The second-order TVD scheme appears to be not stable for the 
chosen limiter and entropy fix. Increasing the order of the WENO scheme from fifth to 
seventh-order gives slightly better results. In addition, the regular WENO-LLF schemes 
are slightly more accurate than their positive counterparts. WEN05-Roe performs with 
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Figure 1. Density plot for the Noh-ID problem obtained by UPWIND (line 1), WEN05-LLF 
(line 2), WEN05-Roe (line 3), WEN05GP (line 4), WEN05P (line 5), WEN07-LLF (line 6), 
WEN07P (line 7) on a grid N — 267. Reference solution: line 8. The right figure is zoomed in 
the vicinity of x = 0. 



Figure 2. Density plot for the Noh-ID problem obtained by WEN05-Roe on a grid N = 267 
with different values of the entropy fix parameter: e = 0.20 (line 1), e = 0.25 (line 2), e = 0.30 
(line 3), e = 0.35 (line 4). Reference solution: line 5. 


the same accuracy as WEN05-LLF using a large entropy fix. Its accuracy can be im- 
proved by using the appropriate entropy fix parameter e. Figure 2 shows the results by 
WEN05-Roe with different values of e. The value e = 0.3 produces an error close to zero 
in the vicinity of x = 0 instead of an error of 0.5% obtained by the scheme with a value 
of e = 0.25. 

The 2D slices at z = 0 obtained for the 3D case by UPWIND, WEN05P, WEN07P, 
WEN05PH, WEN07PH and WEN05PH-Roe are shown in Figure 3. The grid size 
is 134 x 134 x 134. Note that regular WEN05-LLF and WEN05-Roe are not stable 
for the Noh 3D problem due to computed negative pressure/density. The results by 
WEN05P and WEN07P also contain some points with small negative pressure values. 
This can be fixed by applying an adaptive time-step that uses the positivity-preservation 
condition. The results obtained using WEN05PH, WEN05PH-Roe and WEN07PH with 
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Figure 3. Density contours for the Noh-3D problem on a grid 134 x 134 x 134, slice z = 0. 
Top row: UPWIND, WEN05P and WEN07P. Bottom row: WEN05PH-Roe, WEN05PH and 
WEN07PH. 



Figure 4. Comparison of the results for the Noh-3D problem obtained by UPWIND (line 1), 
WEN05P (line 2), WEN07P (line 3), WEN05PH (line 4) and WEN05PH-Roe (line 5) with 
the reference solution (line 6) for the slice y = 0, z — 0 on a grid 134 x 134 x 134. 
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Figure 5. Density contours for the Sedov-3D problem on a grid 128 x 128 X 128 (slice z = 0). 
Top row: WEN05-LLF, WEN05Ph+split and WEN05PH. Bottom row: WEN05P, WEN07P 
and WEN05PH-Roe. 

the same CFL exhibit positive density and pressure values. In the case of computation 
using WEN05PH-Roe, a high value e = 0.5 has been used to achieve an acceptable 
performance. Figure 4 shows a comparison of the results obtained using the same methods 
with the reference solution for the slice y = 0, z = 0. The increase of the scheme order 
from WEN05P to WEN07P does not improve the quality of the results for this problem. 
This is an expected behavior for cases where solutions are almost constant apart from 
discontinuities. Using Roe flux for WEN05PH leads to significant oscillations. 

3.2. The Sedov problem 

The second test case is the 3D spherical Sedov blast wave (Sedov 1959). The initial 
conditions are p = 1, u = 0 and e = 0.1528415451 exp(— i7 2 /i7o)/i7jj, where R = 
^ x 1 + y 2 + z 2 and Rq = 2/ h. The grid size is h = 0.02 and 7 = 1.4. The density and tem- 
perature contours are shown in Figs. 5 and 6 . The deviation from spherical symmetry for 
positive schemes is bigger than for regular WENO, which is well observed on the temper- 
ature contours. However, the solution obtained using standard WENO schemes contains 
some points with small negative pressure values, whereas, WENOP and WENOPfi ob- 
tain all positive pressure and density values. Switching to Roe’s flux causes even further 
deviation from symmetry (see results obtained using WEN05PH-Roe). Note that the 
density contours obtained by all of these methods look very similar. 

3.3. Mach 2000 jet 

The same Mach 2000 jet problem as that in Zhang & Shu (2010) is considered here with 
7 5== 5/3. The computational domain is [0,1] x [—0.25,0.25]. The initial flow condition 
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Figure 6. Temperature contours for the Sedov-3D problem on a grid 128 x 128 X 128 (slice 
z = 0). Top row: WEN05-LLF, WEN05Pfi+split and WEN05PH. Bottom row: WEN05P, 
WEN07P and WEN05PH-Roe. 


is the ambient gas with (p,u,v,p) = (0.5,0,0,0.4127). The boundary conditions for the 
right, top and bottom are outflow. For the left boundary, (p,u,v,p) = (5,800,0,0.4127) 
if y G [—0.05,0.05] and (p,u,v,p) = (0.5,0,0,0.4127) otherwise. The terminal time is 
0.001. The speed of the jet is 800, which is around Mach 2100 with respect to the sound 
speed in the jet gas. 

For this problem a very small initial CFL value is required for high-order computation 
(about 0.01). For this reason, a variable time-step control is used in the computation. 
After each RK stage the solution is tested using the positivity condition. If the condition 
is not satisfied, the time-step is divided by a factor of 2 and the current RK step is 
repeated again. In this way the computation can be carried out with an average CFL 
value 4 — 8 times larger than the fixed CFL. 

The results by different schemes on the uniform grid 800 x 400 are shown in Figure 
7. Note that regular WEN05-LLF, WEN05-Roe and their nonlinear filter counterparts 
WEN05fi+split exhibit negative pressure. In general, the quality of the results by positive 
WENO schemes increases as the order increases. The obtained results using WEN05PH- 
Roe are a little less dissipative than in the case using WEN05PH-LLF or WEN05P. The 
entropy fix parameter used in WEN05PH-Roe for this case was e = 0.3. 

For the above 3D Noh and the Mach 2000 jet test cases, the nonlinear filter counterparts 
of WEN05P, WEN07P WEN05PH and WEN07PH are not stable due to the wide 
stencil of the wavelet flow sensor. Research is underway to develop a flow sensor with 
a more narrow stencil for such flows. Note that these filter counterparts are stable and 
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Figure 7. From top to bottom: density contours for the Mach 2000 jet problem obtained by 
TVD, WEN05GP, WEN05P, WEN05PH-Roe and WEN07P on a grid 800 x 400. Scales are 
logarithmic. 
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p 1.10546 kg/m 3 
T 6000 K 
p 12.7116 MPa 


p 3.0964 x 10 ~ 4 kg/m 3 
T 300 K 
p 26.771 Pa 


Y He 0.9856 
Y N2 0.0144 


Yq 2 0.21 
Y N2 0.79 


Table 1. High- (left) and low- (right) pressure region initial data 


exhibit better accuracy than its non-filter counterparts for problems that do not have 
the aforementioned special extreme condition. 

3.4. 13 species ID EAST simulation 

The computational domain has a total length of 8.5m. The left part of the domain with 
length 0.1m is a high-pressure region. The right part of the domain with length 8.4m is 
a low-pressure region. The gas mixture consists of 13 species: 


The initial conditions of the high- and low-pressure regions are listed in the Table 1. For 
the left-side boundary the Euler (slip) wall condition is applied, and for the right-side, 
the zero gradient condition is applied for all variables. As mentioned before, EAST is 
a high-temperature and high Mach number viscous non-equilibrium flow consisting of 
13 species. In addition, as most common shock- capturing schemes have been developed 
for problems without source terms, when applied to problems with nonlinear and/or 
stiff source terms these methods can result in spurious solutions, even when solving a 
conservative system of equations with a conservative scheme. Here only selected schemes 
are used for this simulation. 

Figure 8 shows the results from the computation using the Harten-Yee second-order 
TVD scheme (Yee 1989; Yee et al. 1990) for four grids with Ax = 10 -3 m, 5x 10 -4 m, 5x 
10 -5 m and 2.5 x 10 -5 m at time t en d = 0.325 x 10 -4 sec. One can observe a significant 
shift in the shear (left discontinuity) and the shock (right discontinuity) locations as the 
grid is refined. The distance between the shear and the shock shrinks as the grid is refined. 
The difference between shock locations obtained on the grids with Ax = 5 x 10 -5 m and 
2.5 x 10 -5 m is less than 0.3%. Thus the solution using Ax = 5 x 10 -5 m can be considered 
as the reference solution. 

The left subfigure of Figure 9 shows a comparison among five methods obtained on a 
coarse grid ( Ax = 10 -3 m) with the reference solution. The scheme’s labels are defined 
as follows: 

• ACMTVDfi: Second-order central base scheme using ACM flow sensor. See Yee et al. 
(1999) for further information on filter schemes. 

• WEN05-llf: Fifth-order WENO (WEN05) using the local Lax-Friedrichs flux. 

• WEN05P-llf: Positive WEN05 of Zhang & Shu (2012) using the local Lax-Friedrichs 
flux. 

• WEN05PH-llf: Positive WEN05 of Hu et al. (2012) using the local Lax-Friedrichs 
flux. 

The right subfigure of Figure 9 shows a comparison of ACMTVDfi using a different 
weight ft parameter of the ACM flow sensor. The smaller the ft, the smaller the amount of 
TVD dissipation that is used. Among the considered schemes, Figure 9 indicates that the 
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X, m 

Figure 8. 13 species ID EAST problem: Second-order Harten-Yee TVD simulation for three 
grids: Ax = 10 -3 m (line 1), 5 X 10 -4 m (line 2), 5 x 10 -5 m (line 3), 2.5 x 10 -5 m (line 4), 
and Tend — 0.325 X 10 _4 s, with CFL — 0.8. 

least dissipative scheme predicts the shear and shock locations best when compared with 
the reference solution. The results indicate that ACMTVDfi is slightly more accurate than 
WEN05-llf. This is due to the fact that ACMTVDfi reduces the amount of numerical 
dissipation away from high gradient regions. Using the subcell resolution method of 
Wang et al. (2012) for one reaction case by applying it to only one of the reactions 
in this multireaction flow does not improve the performance over standard schemes. 
Further research on the generalization of subcell resolution to multi-reactions needs to 
be explored. 


4. Summary 

The positivity-preserving schemes produce more stable behavior than regular WENO 
for the considered problems. The scheme by Hu et al (2012) achieves slightly better 
positivity preservation than the scheme by Zhang & Shu (2012) using the same CFL 
number. These positivity-preserving schemes also are more diffusive than their standard 
WENO counterparts. Accuracy can be improved when using Roe’s flux with the Hu et 
al. scheme instead of using Lax-Friedrichs, which is required by the Zhang-Shu scheme. 
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Figure 9. ID, 13 species EAST problem: Comparison among methods using 601 point grids 
with CFL = 0.6 and t en d — 3.25 x 10 -5 sec. Left subfigure: Reference solution (TVD on a 10, 001 
point grid) (line 1), TVD (line 2), ACMTVDh (TVDfi) using k = 0.5 (line 3), WEN05-llf (line 
4), WEN05P-llf (line 5), WENOPH-llf (line 6). Right: ACMTVDh, k = 0.15 (line 2), k = 0.2 
(line 3), k = 0.3 (line 4), a = 0.5 (line 5), n = 1 (line 6). See text for method notation. 


Department of Energy by Lawrence Livermore National Laboratory under Contract DE- 
AC52-07NA27344. 


REFERENCES 

Griffiths, D., Stuart, A. & Yee, H. 1992 Numerical wave propagation in hyperbolic 
problems with nonlinear source terms. SIAM J. Numer. Anal. 29 , 1244-1260. 

Hu, X. Y., Adams, N. A. & Shu, C.-W. 2012 Positivity-preserving flux limiters for 
high-order conservative schemes. arXiv:1203.1540v4, preprint submitted to Elssvier. 

Jiang, G. S. & Shu, C. W. 1996 Efficient implementation of weighted ENO schemes. 
J. Comp. Phys. 126 , 202-228. 

Lafon, A. & Yee, H. 1996 Dynamical approach study of spurious steady-state numer- 
ical solutions for nonlinear differential equations, part III: The effects of nonlinear 
source terms in reaction-convection equations. Comput. Fluid Dyn. 6 , 1-36. 

LeVeque, R. & Yee, H. C. 1990 A study of numerical methods for hyperbolic conser- 
vation laws with stiff source terms. J. Comp. Phys. 86 , 187-210. 

Noh, W. F. 1987 Errors for calculations of strong shocks using artificial viscosity and 
an artificial heat flux. J. Comput. Phys. 72 , 78-120. 

Sedov, L. I. 1959 Similarity and Dimensional Methods in Mechanics. New York: Aca- 
demic Press. 

Shu, C. W. & Osher, S. 1988 Efficient implementation of essentially nonoscillatory 
shock-capturing schemes. J. Comp. Phys. 77 , 439471. 

Wang, W., Shu, C., Yee, H. C. & Sjogreen, B. 2012 High order finite difference 
methods with subcell resolution for advection equations with stiff source terms. J. 
Comput. Phys. 231 , 190-214. 

Yee, H. C. 1989 A class of high-resolution explicit and implicit shock- capturing methods. 
VKI lecture series 1989-04. 

Yee, H. C., Klopfer, G. H. & Montagne, J.-L. 1990 High-resolution shock- 


12 


D.V. Kotov, H.C. Fee AND B. Sjogreen 

capturing schemes for inviscid and viscous hypersonic flows. J. Comput. Phys. 88 , 
31-61. 

Yee, H. C., Kotov, D. V., Wang, W. & Shu, C.-W. 2012 Spurious behavior of 
shock-capturing methods: Problems containing stiff source terms and discontinuities. 
In Proceedings of the ICCFD7. The Big Island, Hawaii. 

Yee, H. C., Sandham, N. & Djomehri, M. 1999 Low dissipative high order shock- 
capturing methods using characteristic-based filters. J. Comput. Phys. 150 , 199-238. 

Yee, H. C. & Sjogreen, B. 2007 Development of low dissipative high order filter 
schemes for multiscale Navier-Stokes/MHD systems. J. Comput. Phys. 225 , 910- 
934. 

Yee, H. C. & Sjogreen, B. 2010 High order filter methods for wide range of com- 
pressible flow speeds. In Proc. of ASTRONUM-2010 . San Diego, Calif, expanded 
version submitted to Computers & Fluids. 

Zhang, X. & Shu, C.-W. 2010 On positivity preserving high order discontinuous 
galerkin schemes for compressible Euler equations on rectangular meshes. J. Com- 
put. Phys. 229 , 8918-8934. 

Zhang, X. & Shu, C.-W. 2012 Positivity-preserving high order finite difference WENO 
schemes for compressible Euler equations. J. Comput. Phys. 231 , 2245-2258. 


